#ifndef UTILSPART_H_INCLUDED
#define UTILSPART_H_INCLUDED

void setBox(particula *part, parametros *param)
{
    param->vol=pi*param->npart/param->phi/6.;
    #ifdef CONFINED
        if (param->h0>0.){
            ///La altura de la caja es tamayo caja menos 2a
            param->lbox[0]=pow(param->vol/param->h0,1./3.);
            param->lbox[2]=param->h0*param->lbox[0];
            param->lbox[1]=param->lbox[0];
        }
        else{
            param->lbox[0]=pow(param->vol,1./3.);
            param->lbox[1]=param->lbox[2]=param->lbox[0];
        }
    #else
        param->lbox[0]=pow(param->vol,1./3.);
        param->lbox[1]=param->lbox[2]=param->lbox[0];
    #endif
}

void boundary(particula *part, parametros param)
{
	if(part->pos[0]>param.lbox[0]) part->pos[0]-=param.lbox[0];
	if(part->pos[0]<0.) part->pos[0]+=param.lbox[0];

	if(part->pos[1]>param.lbox[1]) part->pos[1]-=param.lbox[1];
	if(part->pos[1]<0.) part->pos[1]+=param.lbox[1];

	#ifndef CONFINED
        if(part->pos[2]>param.lbox[2]) part->pos[2]-=param.lbox[2];
        if(part->pos[2]<0.) part->pos[2]+=param.lbox[2];
	#endif
}

void printBinary(particula *part, parametros param, double sim_time, int step)
{
    double fx,fy,fz,vx,vy,vz;
    FILE *fe;
    /// Abrir el fichero de entrada en lectura y binario
    fe = fopen("config.bin", "ab");
    fwrite(&sim_time,sizeof(double),1,fe);
    fwrite(&step,sizeof(int),1,fe);
    fwrite(&param.npart,sizeof(int),1,fe);
    fwrite(&param.lbox[0],sizeof(double),1,fe);
    fwrite(&param.lbox[1],sizeof(double),1,fe);
    fwrite(&param.lbox[2],sizeof(double),1,fe);
    for (int i=0;i<param.npart;i++)
    {
        ///Imprimimos las posiciones, velocidades, fuerzas y diametros de todas las particulas
        vx=part[i].vel[0];
        vx=part[i].vel[0];
        vx=part[i].vel[0];

        fx=part[i].force[0];
        fx=part[i].force[0];
        fx=part[i].force[0];

        for (int j=0;j<3;j++)
            fwrite(&part[i].pos[j],sizeof(double),1,fe);
        fwrite(&vx,sizeof(double),1,fe);
        fwrite(&vy,sizeof(double),1,fe);
        fwrite(&vz,sizeof(double),1,fe);
        fwrite(&fx,sizeof(double),1,fe);
        fwrite(&fy,sizeof(double),1,fe);
        fwrite(&fz,sizeof(double),1,fe);
        fwrite(&part[i].diametro,sizeof(double),1,fe);
    }

    fclose(fe);
}

void loadConfig(particula *part, parametros param)
{
    void boundary(particula *part, parametros param);
	ifstream MiArchivo ("system_ini.dat",ios::in);
	char line[150];
	double dummy2;
	double rr;
	if(MiArchivo.fail())
	{
		cout<<"No se pudo encontrar el archivo de reinicio "<<endl;
		exit(1);
	}
	if (MiArchivo.is_open())
	{
		///Importamos las posiciones y velocidades
		for(int i=0;i<param.npart;i++) {
			for(int j=0;j<3;j++){
				MiArchivo>>rr;
				part[i].pos[j]=rr;
				part[i].posOld[j]=rr;
			}
			for(int j=0;j<3;j++){
			    MiArchivo>>rr;
				part[i].vel[j]=rr;
			}
			MiArchivo>>rr;
            part[i].diametro=rr;

			boundary(&part[i], param);
			part[i].mass=pow(part[i].diametro,3.0);
			part[i].gamma=0.;
		}
		MiArchivo.close();
	}
}

void printConfigBackup(parametros param, particula *part)
{
	ofstream MiArchivo ("system_bk.dat",ios::out);

	if (MiArchivo.is_open()) {
        for(int i=0;i<param.npart;i++) {
            MiArchivo<<scientific<<part[i].pos[0]<<setw(20)<<part[i].pos[1]<<setw(20)<<part[i].pos[2]<<
            setw(20)<<part[i].vel[0]<<setw(20)<<part[i].vel[1]<<setw(20)<<part[i].vel[2]<<setw(20)<<part[i].diametro<<endl;
        }
	}
}

void SetRandomPart( particula *part, parametros *param)
{
    double distancia2(particula part1, particula part2, parametros *param);
    int nok=0;
    bool good;
    double r2;
    double dmin;
    double Vpart=0.0;
    double u1,u2;
    double daux;
    for(int i=0;i< param->npart;i++)
    {
        u1=double(rand())/RAND_MAX;
        u2=double(rand())/RAND_MAX;
        part[i].diametro=exp(param->poly*sqrt(-2.*log(u1))*cos(2.*pi*u2));
        Vpart+=1./6.*pi*pow(part[i].diametro,3.0);
    }
    param->vol=Vpart/param->phi;

    param->lbox[0]=pow(param->vol/param->h0,1./3.);
    param->lbox[1]=param->lbox[0];
    param->lbox[2]=param->h0*param->lbox[0];
    dmin=param->lbox[0];

	for(int i=0;i< param->npart;i++)
	{
        part[i].pos[0]=param->lbox[0]*(double(rand())/RAND_MAX);
        part[i].pos[1]=param->lbox[1]*(double(rand())/RAND_MAX);
        part[i].pos[2]=param->lbox[2]*(double(rand())/RAND_MAX);

		///Comprobamos solapamiento con otras particulas
		good=true;
        daux=1.e7;
		for(int j=0;j<nok;j++){
		    ///Recorremos todas las particulas que ya esten metidas en la lista
		    if(good){
		        r2=distancia2(part[i],part[j],param);
                if(r2<1.1*0.5*(part[i].diametro+part[j].diametro)){
                    ///Hay solapamiento entre las particulas i y j. No aceptamos la particula i
                    good=false;
                    i--;
                }
                else
                {
                    if(r2<daux) daux=r2;
                }
		    }
		}
		if(good){

            if( dmin>daux) dmin=daux;
            part[i].mass=part[i].diametro*part[i].diametro*part[i].diametro;
            for(int j=0;j<3;j++) part[i].posOld[j]=part[i].pos[j];
            part[i].gamma=0.;
            ///Si pasa por aqui la particula no se solapa con ninguna otra --> se acepta
            nok++;
		}
	}
	cout<<"Particulas colocadas en la caja: (dmin: "<<dmin<<")"<<endl;
}

void setVelocidades(particula *part, parametros param)
{
    void rescaleVel(particula *part, parametros param);
    double f,v2=0;
    for(int i=0;i<param.npart;i++)
    {
        v2=0.0;
        for(int j=0;j<3;j++)
        {
            part[i].vel[j]=3.0*(double(rand())/RAND_MAX-0.5);
            v2+=part[i].vel[j]*part[i].vel[j];
        }
    }
    rescaleVel(part, param);
}

void rescaleVel(particula *part, parametros param)
{
    double v2, f;
    for(int i=0;i<param.npart;i++)
    {
        v2=0.0;
        for(int j=0;j<3;j++)
        {
            v2+=part[i].vel[j]*part[i].vel[j];
        }
        f=3./v2;
        for(int j=0;j<3;j++)
        {
            part[i].vel[j]=part[i].vel[j]*sqrt(f);
        }
    }
}

double distancia2(particula part2, particula part1, parametros *param)
{
    double dis2=0;
    int dim;
    double ddr[3];
    int i;

    for (i=0;i<2;i++)
    {
        ddr[i]=part1.pos[i]-part2.pos[i];
        if(ddr[i]>0.5*param->lbox[i])
            ddr[i]-=param->lbox[i];
        if(ddr[i]<-0.5*param->lbox[i])
            ddr[i]+=param->lbox[i];
        dis2+=ddr[i]*ddr[i];
    }

    ddr[2]=part1.pos[2]-part2.pos[2];
    #ifndef CONFINED
        if(ddr[2]>0.5*param->lbox[2])
            ddr[2]-=param->lbox[2];
        if(ddr[2]<-0.5*param->lbox[2])
            ddr[2]+=param->lbox[2];
        dis2+=ddr[i]*ddr[i];
    #endif

    for (i=0;i<3;i++) param->dr[i]=ddr[i];
    return  sqrt(dis2);
}

#endif // UTILSPART_H_INCLUDED
